pro plot_cpu_usage,cpufile,makefile,gadget2=gadget2,non_cosmo=non_cosmo,max_y=max_y,outfile=outfile,constrain_x=constrain_x

;+
; NAME:
;       PLOT_CPU_USAGE
;
; PURPOSE:
;       This routine reads in cpu.txt files from Gadget2/3 and plots
;       the various components
;
; CATEGORY:
;       
;
; CALLING SEQUENCE:
;       PLOT_CPU_USAGE(FILE,MAKEFILE[,GADGET2=GADGET2,NON_COSMO=NON_COSMO,MAX_Y=MAX_Y])
;
; INPUTS:
;       FILE:      Name of the cpu file to read in.
;       MAKEFILE:  Name of the Makefile used to compile the executable,
;                  which is parsed to retrieve the value of PMGRID
;
; OPTIONAL INPUTS:
;       GADGET2:   Instructs routine to expect the (very different)
;                  cpu.txt format produced by public Gadget2
;       NON_COSMO: By default the time variable is considered to be 
;                  the cosmological expansion factor; as such the 
;                  routine adds an upper axis for redshift. For 
;                  non-cosmological runs, the time variable is
;                  converted from internal units to Myr and the upper
;                  axis is omitted.
;       MAX_Y:     Upper limit of the y-axis. Defaults to 110% of the 
;                  actual max, but it is often useful to specify a
;                  fixed value to enable direct comparison between two
;                  runs.
;
; OUTPUTS:
;       cpu.pdf    Routine produces a .ps file which is subsequently
;                  converted to .pdf by ps2pdf, and the .ps is deleted
;
; RESTRICTIONS:
;
;
; PROCEDURE:
;
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;       Started by Rob Crain, Swinburne, 16-Feb-2010
;-

if(n_params() LT 2)then begin
   print,' usage: plot_cpu_usage,cpufile,makefile[,gadget2=gadget2,non_cosmo=non_cosmo,max_y=max_y]'
   stop
endif

junkdir = '/nfs/cluster/blake/rcrain/junk/'

if keyword_set(gadget2) then begin
   print,' assuming Gadget2 cpu file'
endif else print,' assuming Gadget3-BG cpu file'
if keyword_set(non_cosmo) then begin
   print,' assuming non-cosmological simulation'
endif else print,' assuming cosmological simulation'

; get size of pmgrid from Makefile
cmd = "grep DPMGRID "+makefile
spawn,cmd,pmstatement
; test whether -DPMGRID has been commented
tmp = strpos(pmstatement,'#')
if(tmp EQ -1)then begin
   cmd = "grep DPMGRID "+makefile+" | awk '{print $3}' "
   spawn,cmd,pmstatement   
   pmtrim      = 2
endif else begin
   pmstatement = 'PMGRID unused'
   pmtrim      = 0
endelse

print," -found PMGRID:",pmstatement

; find line number of last occurence of complete time-step
print," -finding last complete time-step..."
cmd = 'grep -n misc '+cpufile+' | tail -1'
spawn,cmd,result
tmp = strsplit(result,':',/extract)
cn  = tmp[0]
cmd = 'head -'+cn+' '+cpufile+' > '+junkdir+'/tmp.file'
spawn,cmd,result
print,"  -found "+cn+" complete time-steps"

; grep out the various elements
cmd = "grep -i Time "+junkdir+"/tmp.file -m 1 | awk '{print $6}' | head -1 "
spawn,cmd,result 
ncpu = long(result)
print,' -ncpu = ',ncpu

print," -extracting components..."

cmd = "grep -i Time "+junkdir+"/tmp.file | awk '{print $4}' > "+junkdir+"/tmp.time ; cat "+junkdir+"/tmp.time | wc -l"
spawn,cmd,result,err 
time = fltarr(long(reform(result)))
print,' time:     found nlines = ',result
;
cmd = "grep -i total "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.total ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
total    = fltarr(long(reform(result)))
allother = fltarr(long(reform(result)))
print,' total:    found nlines = ',result
;
cmd = "grep -i treegrav "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.treegrav ; cat "+junkdir+"/tmp.treegrav | wc -l"
spawn,cmd,result 
treegrav  = fltarr(long(reform(result)))
treeother = fltarr(long(reform(result)))
print,' total:    found nlines = ',result
;
cmd = "grep -i treewalk "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.treewalk ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
treewalk = fltarr(long(reform(result)))
print,' treewalk: found nlines = ',result
;
cmd = "grep -i treecomm "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.treecomm ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
treecomm = fltarr(long(reform(result)))
print,' treecomm: found nlines = ',result
;
cmd = "grep -i sph "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.sph ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
sph      = fltarr(long(reform(result)))
sphother = fltarr(long(reform(result)))
print,' sph:      found nlines = ',result
;
cmd = "grep -i density "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.density ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
density = fltarr(long(reform(result)))
print,' density:  found nlines = ',result
;
cmd = "grep -i denscomm "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.denscomm ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
denscomm = fltarr(long(reform(result)))
print,' denscomm: found nlines = ',result
;
cmd = "grep -i hydcomm "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.hydcomm ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
hydcomm = fltarr(long(reform(result)))
print,' hydcomm:  found nlines = ',result
;
cmd = "grep -i pmgrav "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.pmgrav ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
pmgrav = fltarr(long(reform(result)))
print,' pmgrav:   found nlines = ',result
;
cmd = "grep -i domain "+junkdir+"/tmp.file | awk '{print $2}' > "+junkdir+"/tmp.domain ; cat "+junkdir+"/tmp.total | wc -l"
spawn,cmd,result 
domain = fltarr(long(reform(result)))
print,' domain:   found nlines = ',result

; read each value
close,/all
;
openr,1,junkdir+'/tmp.time'
readf,1,time
close,1
;
openr,1,junkdir+'/tmp.total'
readf,1,total
close,1
;
openr,1,junkdir+'/tmp.treegrav'
readf,1,treegrav
close,1
;
openr,1,junkdir+'/tmp.treewalk'
readf,1,treewalk
close,1
;
openr,1,junkdir+'/tmp.treecomm'
readf,1,treecomm
close,1
;
openr,1,junkdir+'/tmp.sph'
readf,1,sph
close,1
;
openr,1,junkdir+'/tmp.density'
readf,1,density
close,1
;
openr,1,junkdir+'/tmp.denscomm'
readf,1,denscomm
close,1
;
openr,1,junkdir+'/tmp.hydcomm'
readf,1,hydcomm
close,1
;
openr,1,junkdir+'/tmp.pmgrav'
readf,1,pmgrav
close,1
;
openr,1,junkdir+'/tmp.domain'
readf,1,domain
close,1

; if non-cosmological, we need to convert internal units into Myr
if keyword_set(non_cosmo)then time *= 500.0/0.372449

; convert all into hours
total    /= 3.6e3
treegrav /= 3.6e3
treewalk /= 3.6e3
treecomm /= 3.6e3
sph      /= 3.6e3
density  /= 3.6e3
denscomm /= 3.6e3
hydcomm  /= 3.6e3
pmgrav   /= 3.6e3
domain   /= 3.6e3

; generate 'other' values
allother   = total    - treegrav - pmgrav   - sph     - domain
sphother   = sph      - hydcomm  - denscomm - density
treeother  = treegrav - treewalk - treecomm

;; find occasions where snapshot-based restarts have reset the
;; wallclock, and add on the necessary offset.
;ntot = n_elements(total)
;for i=0l,long(ntot)-2l do begin
;    if(total[i+1] LT total[i])then begin
;        print,' -adding offset to total of:',total[i]-total[i+1]
;        total[i+1<ntot-1:ntot-1]    += total[i]-total[i+1]
;        treegrav [i+1<ntot-1:ntot-1] += treegrav[i]-treegrav[i+1]
;        treewalk [i+1<ntot-1:ntot-1] += treewalk[i]-treewalk[i+1]
;        treecomm [i+1<ntot-1:ntot-1] += treecomm[i]-treecomm[i+1]
;        sph      [i+1<ntot-1:ntot-1] += sph[i]-sph[i+1]
;        density  [i+1<ntot-1:ntot-1] += density[i]-density[i+1]
;        denscomm [i+1<ntot-1:ntot-1] += denscomm[i]-denscomm[i+1]
;        hydcomm  [i+1<ntot-1:ntot-1] += hydcomm[i]-hydcomm[i+1]
;        pmgrav   [i+1<ntot-1:ntot-1] += pmgrav[i]-pmgrav[i+1]
;        domain   [i+1<ntot-1:ntot-1] += domain[i]-domain[i+1]
;        allother [i+1<ntot-1:ntot-1] += allother[i]-allother[i+1] 
;        sphother [i+1<ntot-1:ntot-1] += sphother[i]-sphother[i+1]
;        treeother[i+1<ntot-1:ntot-1] += treeother[i]-treeother[i+1]
;    endif
;endfor
;
@set_rcrain_colours

set_plot,'ps'
if keyword_set(outfile) then begin
   myplot = outfile
endif else myplot = 'timing_g3'
myplot_ps  = myplot+'.ps'
myplot_pdf = myplot+'.pdf'
device,/color,filename=myplot_ps,xsize=18.5,ysize=16,encapsulated=1    
!p.font = -1    


ylab = TEXtoIDL("Wallclock [hrs]")
dummy = fltarr(1)

if keyword_set(non_cosmo) then begin
   xstyle = 1
   if keyword_set(constrain_x)then begin
       xlims = minmax(time)
   endif else xlims  = [0,500]
   xlab = TEXtoIDL("Time [Myr]")
endif else begin
   xstyle = 9
   if keyword_set(constrain_x)then begin
       xlims = minmax(time)
   endif else xlims  = [0,1]
   xlab = TEXtoIDL("Expansion Factor")
endelse
if keyword_set(max_y)then begin
   ylims = [0,max_y]
endif else ylims=[0,1.1*max(total)]


plot,time,total,xs=xstyle,ys=1,xr=xlims,yr=ylims,$
     thick=3.0,xthick=3.0,ythick=3.0,xtit=xlab,ytit=ylab,charsize=1.3,charthick=3.0,$
     position=[0.130,0.1,0.95,0.90]

if not(keyword_set(non_cosmo)) then begin
   z_vals_we_want = [0.,1.,2.,5.,10.,75.]
   a_vals_we_want = 1.0 / (z_vals_we_want+1.0)
   x2lab = TEXtoIDL("Redshift")
   axis,xaxis=1,xstyle=1,color=black,xtit=x2lab,$
        xthick=3.0,xticks=5,xtickv=a_vals_we_want,xtickname=blanknames,xtickformat='z_from_a',$
        charsize=1.3,charthick=3.0
endif



; draw and shade these regions
data = total 
xpoly = [time,reverse(time)]
ypoly = [replicate(0,n_elements(time)),reverse(data)]
PolyFill, xpoly, ypoly, Color=orange, noclip=0

data = total - allother
xpoly = [time,reverse(time)]
ypoly = [replicate(0,n_elements(time)),reverse(data)]
PolyFill, xpoly, ypoly, Color=cyan, noclip=0

data = data - domain
xpoly = [time,reverse(time)]
ypoly = [replicate(0,n_elements(time)),reverse(data)]
PolyFill, xpoly, ypoly, Color=red, noclip=0

data = data - sph
xpoly = [time,reverse(time)]
ypoly = [replicate(0,n_elements(time)),reverse(data)]
PolyFill, xpoly, ypoly, Color=green, noclip=0

data = data - pmgrav
xpoly = [time,reverse(time)]
ypoly = [replicate(0,n_elements(time)),reverse(data)]
PolyFill, xpoly, ypoly, Color=blue, noclip=0

data = total 
oplot,time,data,color=black,thick=5.0
data -= allother
oplot,time,data,color=black,thick=5.0
data -= domain
oplot,time,data,color=black,thick=5.0
data -= sphother
oplot,time,data,color=black,thick=1.0,linestyle=1
data -= hydcomm
oplot,time,data,color=black,thick=1.0,linestyle=1
data -= denscomm
oplot,time,data,color=black,thick=1.0,linestyle=1
data -= density
oplot,time,data,color=black,thick=5.0
data -= pmgrav
oplot,time,data,color=black,thick=5.0
data -= treeother
oplot,time,data,color=black,thick=1.0,linestyle=1
data -= treecomm
oplot,time,data,color=black,thick=1.0,linestyle=1
data -= treewalk

; plot legend
if keyword_set(gadget2) then begin
   code_label = 'Gadget2'
endif else code_label = 'Gadget3-BG'
xyouts,0.45,0.820,TEXtoIDL(code_label)                            ,charsize=1.0,charthick=3.0,/normal
xyouts,0.45,0.780,TEXtoIDL('N_{CPU} = ')+strtrim(string(ncpu),2)  ,charsize=1.0,charthick=3.0,/normal
xyouts,0.45,0.740,strtrim(strmid(pmstatement,pmtrim))             ,charsize=1.0,charthick=3.0,/normal
;
xyouts,0.18,0.820,'Total'        ,charsize=1.0,charthick=3.0,/normal
xyouts,0.18,0.780,'Other'        ,charsize=1.0,charthick=3.0,/normal
xyouts,0.18,0.740,'Domain'       ,charsize=1.0,charthick=3.0,/normal
xyouts,0.18,0.700,'SPH'          ,charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.670,'Other'        ,charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.640,'Hydro comms'  ,charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.610,'Density comms',charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.580,'Density calc' ,charsize=1.0,charthick=3.0,/normal
xyouts,0.18,0.540,'PM gravity'   ,charsize=1.0,charthick=3.0,/normal
xyouts,0.18,0.500,'Tree gravity' ,charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.470,'Other'        ,charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.440,'Tree comm'    ,charsize=1.0,charthick=3.0,/normal
xyouts,0.20,0.410,'Tree walk'    ,charsize=1.0,charthick=3.0,/normal

plots,[0.32,0.38],[0.830,0.830],color=black ,thick=15.0,/normal
plots,[0.32,0.38],[0.790,0.790],color=orange,thick=15.0,/normal
plots,[0.32,0.38],[0.750,0.750],color=cyan  ,thick=15.0,/normal
plots,[0.32,0.38],[0.710,0.710],color=red   ,thick=15.0,/normal
plots,[0.32,0.38],[0.550,0.550],color=green ,thick=15.0,/normal
plots,[0.32,0.38],[0.510,0.510],color=blue  ,thick=15.0,/normal

device,/close
print,'ps2pdf -dEPSCrop '+myplot_ps+' '+myplot_pdf
spawn,'ps2pdf -dEPSCrop '+myplot_ps+' '+myplot_pdf
spawn,'\rm '+myplot_ps

cmd = "\rm "+junkdir+"/tmp.*"
spawn,cmd,result

stop
end



function z_from_a,axis, index, value
; given the expansion factor, return redshift
age =  galage(value,1000,H0=73,Omega=0.25,Lambda=0.75,/SILENT)
age = age / 1.0e9

z  = (1.0 / value) - 1.0

if(value eq 0.5)then begin
   zs = '0.5'
endif else begin
   zs = nint(z)
endelse

zs = strtrim(string(zs),2)

RETURN,zs
end
